Nota Bene (NB): This is an old document that is still valid. While I
have added some new data projections, the basic concepts are the same. I hope to clean this up soon.
M.F. 28 May 1997
PROJECTIONS in GrADS
Mike Fiorino
PCMDI/LLNL
fiorino@typhoon.llnl.gov
26 April, 1995
Here's everything you wanted to know about projections in GrADS
but were afraid to ask. I hope this doc is useful; your comments
are appreciated.
Let's get clear about projections in GrADS -- there are TWO
issues: 1) projection of the DATA; and 2) projection of the
DISPLAY.
GrADS supports two types of DATA grids: 1) lon/lat grids (and NOT
necessarily regular, e.g., gaussian); and 2) PREPROJECTED grids.
The latter are DATA ALREADY on a map projection.
GrADS supports four types of PREPROJECTED data: 1) N polar
stereo; 2) S polar stereo; 3) lambert conformal; and 4) (beta)
NMC eta model.
When PREPROJECTED grids are opened in GrADS, bilinear
interpolation constants are calculated and all data are displyed
on an INTERNAL GrADS LAT/LON GRID DEFINED BY the XDEF and YDEF
card in the data description or ".ctl" file (that's why it takes
longer to "open" a preprojected grid data set).
It is very important to point out that the INTERNAL GrADS grid
can be ANY GRID as it is COMPLETELY INDEPENDENT OF THE
PREPROJECTED DATA GRID. Thus, there is nothing stopping you from
displaying preprojected data on a very high res LON/LAT grid
(again, defined in the .ctl by xdef and ydef). In fact, you
could create and open MULTIPLE .ctl files with DIFFERENT
resolutions and/or regions which pointed to the SAME PREPROJECTED
data file.
When you do a "display" (i.e., get a grid of data), the
PREPROJECTED data are bilinearly interpolated to the GrADS
INTERNAL lat/lon grid. For PREPROJECTED scalar fields (e.g., 500
mb heights), the display is adequate and the PRECISION of the
interpolation can be controlled by xdef and ydef to define a
higher spatial resolution grid.
The big virtue of this approach is that all built in GrADS
analytic functions (e.g., aave, hcurl...) continue to work even
though the data were NOT originally on a lon/lat grid. The
downside is that you are not looking directly at your data ON A
GEOGRAPHIC MAP. However, one could always define a .ctl file
which simply opened the data file as i,j data and displayed
WITHOUT the map (set mpdraw off). So, in my opinion, this
compromise is not that limiting even if as a modeler you wanted
to look at the grid -- you just don't get the map background.
PREPROJECTED vector fields are a little trickier, depending on
whether the vector is defined relative to the DATA grid or
relative to the EARTH. For example, NMC polar stereo grids use
winds relative to the DATA grid and thus must be ROTATED to the
INTERNAL GrADS lat/lon grid (again defined in the .ctl file by
the xdef and ydef cards).
The only potential problem with working with PREPROJECTED data
(e.g., lambert conformal model data) is DEFINING the projection
for GrADS. This is accomplished using a pdef card in the data
descriptor ".ctl" file.
---------- POLAR STEREO PREPROJECTED DATA ---------
PREPROJECTED data on a polar stereo projection (N and S) is
defined a la NMC. For the NMC NGM model GRIB data distributed
via anon ftp to nic.fb4.noaa.gov, the pdef card is:
* pdef isize jsize projtype ipole jpole lonref gridinc
pdef 53 45 nps 27 49 -105 190.5
(NOTE: the * in the first column of the .ctl file means a comment...)
where,
ipole and jpole are the (i,j) of the pole and gridinc is the dx in KM!!!
The relevent GrADS source is:
void w3fb04 (float alat, float along, float xmeshl, float orient,
float *xi, float *xj) {
/*
C
C SUBPROGRAM: W3FB04 LATITUDE, LONGITUDE TO GRID COORDINATES
C AUTHOR: MCDONELL,J. ORG: W345 DATE: 90-06-04
C
C ABSTRACT: CONVERTS THE COORDINATES OF A LOCATION ON EARTH FROM THE
C NATURAL COORDINATE SYSTEM OF LATITUDE/LONGITUDE TO THE GRID (I,J)
C COORDINATE SYSTEM OVERLAID ON A POLAR STEREOGRAPHIC MAP PRO-
C JECTION TRUE AT 60 DEGREES N OR S LATITUDE. W3FB04 IS THE REVERSE
C OF W3FB05.
C
C PROGRAM HISTORY LOG:
C 77-05-01 J. MCDONELL
C 89-01-10 R.E.JONES CONVERT TO MICROSOFT FORTRAN 4.1
C 90-06-04 R.E.JONES CONVERT TO SUN FORTRAN 1.3
C 93-01-26 B. Doty converted to C
C
C USAGE: CALL W3FB04 (ALAT, ALONG, XMESHL, ORIENT, XI, XJ)
C
C INPUT VARIABLES:
C NAMES INTERFACE DESCRIPTION OF VARIABLES AND TYPES
C ------ --------- -----------------------------------------------
C ALAT ARG LIST LATITUDE IN DEGREES (<0 IF SH)
C ALONG ARG LIST WEST LONGITUDE IN DEGREES
C XMESHL ARG LIST MESH LENGTH OF GRID IN KM AT 60 DEG LAT(<0 IF SH)
C (190.5 LFM GRID, 381.0 NH PE GRID,-381.0 SH PE GRID)
C ORIENT ARG LIST ORIENTATION WEST LONGITUDE OF THE GRID
C (105.0 LFM GRID, 80.0 NH PE GRID, 260.0 SH PE GRID)
C
C OUTPUT VARIABLES:
C NAMES INTERFACE DESCRIPTION OF VARIABLES AND TYPES
C ------ --------- -----------------------------------------------
C XI ARG LIST I OF THE POINT RELATIVE TO NORTH OR SOUTH POLE
C XJ ARG LIST J OF THE POINT RELATIVE TO NORTH OR SOUTH POLE
C
C SUBPROGRAMS CALLED:
C NAMES LIBRARY
C ------------------------------------------------------- --------
C COS SIN SYSLIB
C
C REMARKS: ALL PARAMETERS IN THE CALLING STATEMENT MUST BE
C REAL. THE RANGE OF ALLOWABLE LATITUDES IS FROM A POLE TO
C 30 DEGREES INTO THE OPPOSITE HEMISPHERE.
C THE GRID USED IN THIS SUBROUTINE HAS ITS ORIGIN (I=0,J=0)
C AT THE POLE IN EITHER HEMISPHERE, SO IF THE USER'S GRID HAS ITS
C ORIGIN AT A POINT OTHER THAN THE POLE, A TRANSLATION IS NEEDED
C TO GET I AND J. THE GRIDLINES OF I=CONSTANT ARE PARALLEL TO A
C LONGITUDE DESIGNATED BY THE USER. THE EARTH'S RADIUS IS TAKEN
C TO BE 6371.2 KM.
C
C ATTRIBUTES:
C LANGUAGE: SUN FORTRAN 1.4
C MACHINE: SUN SPARCSTATION 1+
C*/
static float radpd = 0.01745329;
static float earthr = 6371.2;
float re,xlat,wlong,r;
re = (earthr * 1.86603) / xmeshl;
xlat = alat * radpd;
if (xmeshl>0.0) {
wlong = (along + 180.0 - orient) * radpd;
r = (re * cos(xlat)) / (1.0 + sin(xlat));
*xi = r * sin(wlong);
*xj = r * cos(wlong);
} else {
re = -re;
xlat = -xlat;
wlong = (along - orient) * radpd;
r = (re * cos(xlat)) / (1.0+ sin(xlat));
*xi = r * sin(wlong);
*xj = -r * cos(wlong);
}
}
---------- LAMBERT CONFORMAL PREPROJECTED DATA ---------
The lambert conformal projection (lcc) was set up for the
U.S. Navy's limited area model NORAPS. Thus, to work with YOUR
lcc data you must express YOUR grid in the context of the NAVY
lcc grid. NMC has been able to do this for their AIWIPS grids
and the Navy definition should be general enough for others.
A typical NORAPS lambert-conformal grid is described below,
including the C code which sets up the internal interpolation.
The .ctl file is:
dset ^temp.grd
title NORAPS DATA TEST
undef 1e20
pdef 103 69 lcc 30 -88 51.5 34.5 20 40 -88 90000 90000
xdef 180 linear -180 1.0
ydef 100 linear -10 1.0
zdef 16 levels 1000 925 850 700 500 400 300 250 200 150 100 70 50 30 20 10
tdef 1 linear 00z1jan94 12hr
vars 1
t 16 0 temp
endvars
where,
103 = #pts in x
69 = #pts in y
lcc = lambert-conformal
30 = lat of a ref point
-88 = lon of ref point (E is positive in GrADS, W is negative)
51.5 = i of ref point
34.5 = j of ref point
20 = S true lat
40 = N true lat
-88 = standard lon
90000 = dx in M
90000 = dy in M
GrADS source which maps lon/lat of the GrADS INTERNAL lon/lat
grid to i,j of the PREPROJECTED grid....
/* Lambert conformal conversion */
void ll2lc (float *vals, float grdlat, float grdlon,
float *grdi, float *grdj) {
/* Subroutine to convert from lat-lon to lambert conformal i,j.
Provided by NRL Monterey; converted to C 6/15/94.
c SUBROUTINE: ll2lc
c
c PURPOSE: To compute i- and j-coordinates of a specified
c grid given the latitude and longitude points.
c All latitudes in this routine start
c with -90.0 at the south pole and increase
c northward to +90.0 at the north pole. The
c longitudes start with 0.0 at the Greenwich
c meridian and increase to the east, so that
c 90.0 refers to 90.0E, 180.0 is the inter-
c national dateline and 270.0 is 90.0W.
c
c INPUT VARIABLES:
c
c vals+0 reflat: latitude at reference point (iref,jref)
c
c vals+1 reflon: longitude at reference point (iref,jref)
c
c vals+2 iref: i-coordinate value of reference point
c
c vals+3 jref: j-coordinate value of reference point
c
c vals+4 stdlt1: standard latitude of grid
c
c vals+5 stdlt2: second standard latitude of grid (only required
c if igrid = 2, lambert conformal)
c
c vals+6 stdlon: standard longitude of grid (longitude that
c points to the north)
c
c vals+7 delx: grid spacing of grid in x-direction
c for igrid = 1,2,3 or 4, delx must be in meters
c for igrid = 5, delx must be in degrees
c
c vals+8 dely: grid spacing (in meters) of grid in y-direction
c for igrid = 1,2,3 or 4, delx must be in meters
c for igrid = 5, dely must be in degrees
c
c grdlat: latitude of point (grdi,grdj)
c
c grdlon: longitude of point (grdi,grdj)
c
c grdi: i-coordinate(s) that this routine will generate
c information for
c
c grdj: j-coordinate(s) that this routine will generate
c information for
c
*/
float pi, pi2, pi4, d2r, r2d, radius, omega4;
float gcon,ogcon,ahem,deg,cn1,cn2,cn3,cn4,rih,xih,yih,rrih,check;
float alnfix,alon,x,y;
pi = 4.0*atan(1.0);
pi2 = pi/2.0;
pi4 = pi/4.0;
d2r = pi/180.0;
r2d = 180.0/pi;
radius = 6371229.0;
omega4 = 4.0*pi/86400.0;
/*mf -------------- mf*/
/*case where standard lats are the same */
if(*(vals+4) == *(vals+5)) {
gcon = sin(*(vals+4)*d2r);
} else {
gcon = (log(sin((90.0-*(vals+4))*d2r))
-log(sin((90.0-*(vals+5))*d2r)))
/(log(tan((90.0-*(vals+4))*0.5*d2r))
-log(tan((90.0-*(vals+5))*0.5*d2r)));
}
/*mf -------------- mf*/
ogcon = 1.0/gcon;
ahem = fabs(*(vals+4))/(*(vals+4));
deg = (90.0-fabs(*(vals+4)))*d2r;
cn1 = sin(deg);
cn2 = radius*cn1*ogcon;
deg = deg*0.5;
cn3 = tan(deg);
deg = (90.0-fabs(*vals))*0.5*d2r;
cn4 = tan(deg);
rih = cn2*pow((cn4/cn3),gcon);
deg = (*(vals+1)-*(vals+6))*d2r*gcon;
xih = rih*sin(deg);
yih = -rih*cos(deg)*ahem;
deg = (90.0-grdlat*ahem)*0.5*d2r;
cn4 = tan(deg);
rrih = cn2*pow((cn4/cn3),gcon);
check = 180.0-*(vals+6);
alnfix = *(vals+6)+check;
alon = grdlon+check;
while (alon<0.0) alon = alon+360.0;
while (alon>360.0) alon = alon-360.0;
deg = (alon-alnfix)*gcon*d2r;
x = rrih*sin(deg);
y = -rrih*cos(deg)*ahem;
*grdi = *(vals+2)+(x-xih)/(*(vals+7));
*grdj = *(vals+3)+(y-yih)/(*(vals+8));
}
There are a few GOTCHAs with using PREPROJECTED data:
1) the units in the variable definition for the u and v
components MUST BE 33 and 34 (the GRIB standard) respectively,
e.g.,
u 15 33 u component of the wind at 15 pressure levels
v 15 34 v component of the wind at 15 pressure levels
2) wind rotation is handled for POLAR STEREO (N and S)
PREPROJECTED data, but NOT FOR LAMBERT CONFORMAL, as the Navy
rotates the winds to earth relative. This will have to be added
later......
3) the eta projection is still experimental...
---------- GrADS DISPLAY PROJECTIONS ---------
Now that you hopefully understand GrADS DATA GRIDS, it is time to
discuss DISPLAY PROJECTIONS. Graphics in GrADS are calculated
relative to the INTERNAL GrADS DATA GRID I,J space, TRANSFORMED
to the display device coordinates (e.g., the screen) and then
displayed. That is, the i,j of the graphic element is converted
to lat/lon and then to x,y on the screen via A MAP PROJECTION.
GrADS currently supports four DISPLAY PROJECTIONS: 1) lat/lon (or
spherical); 2) N polar stereo (set mproj nps); 3) S polar stereo;
and 4) the Robinson projection (set lon -180 180, set lat -90 90,
set mproj robinson). As you can probably appreciate, the
i,j-to-lon/lat-to-screen x,y for LON/LAT displays is very simple
and is considerably more complicated for N and S polar stereo
projections.
IN PRINCIPLE, a Lambert Conformal display projection could be
implemented. It just takes work and a simple user interface for
setting up that DISPLAY projection. Actually, the user interface
(i.e., "set" calls) is the most difficult problem...
---------- SUMMARY and PLANS ---------
GrADS handles map projections in two different ways. The first
is PREPROJECTED data where the fields are ALREADY on a projection
(e.g., lambert conformal). It is fairly straightforward to
implement OTHER PREPROJECTED DATA projections and I will be fully
implementing the NMC eta grid both staggered and unstaggered,
"thinned" gaussian grids and the CSU RAMS oblique polar stereo
projection. The second is in how i,j graphics (calculated in
"grid" space) are displayed on a map background. Currently, only
a few basic projections (lon/lat, polar stereo and robinson) are
supported, but perhaps the development group will tackle this
problem. I'm personally quite satisfied with the current display
projections, but then I work with global data sets.
I hope this makes sense; your comments are invited and
encouraged.